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ABSTRACT 

In this paper, we examine some splitting techniques for low Mach number 
Euler flows. We point out shortcomings of some of the proposed methods and 
suggest an explanation for their inadequacy. We then present a symmetric 
splitting for both the Euler and Navier-Stokes equations which removes the 
stiffness of these equations when the Mach number is small. The splitting is 
shown to be stable. 
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INTRODUCTION 


For many computational problems in low speed fluid-dynamics, it has been 
customary to use the incompressible Euler or Navier-Stokes equations. There 
are essentially two reasons for doing this: there is one less variable, since 
the density remains constant, and the stability limit is independent of the 
sound speed. Recently, however, there has been increased interest in studying 
compressibility effects even for low Mach number fluid flows. The compress- 
ible equations, unfortunately, have stiff coefficients due to the disparity in 
the magnitude of the flow velocity and the speed of sound. To overcome this 
difficulty various splitting methods have been proposed to remove the stiff- 
ness from the matrix coefficients of the equations, [3, 4, 7]. Some of these 
methods, however, have not performed as anticipated; in fact, often, for the 
the stipulated stability limits on the time step, the calculations diverged. 

In this paper, we first propose an explanation for this behavior. We 
give examples in the first three sections which show that splittings resulting 
in matrices which are not simultaneously symmetrizable (such as in [7]) may be 
ill-posed at the p.d.e. level. Similar results are presented for some 
explicit numerical schemes, both finite difference and spectral. Thus, the 
intent of these sections is to caution against unrestrained use of splitting 
methods. 

In Section IV, we present a transformation of variables which symmetrizes 
the Euler equations. Under the assumption of low Mach number flow, we are 
able to propose an efficient splitting technique for the compressible equa- 
tions. The resulting algorithm, given both for the Euler and Navier-Stokes" 
equations, is unstiff for the nonlinear field, and the other split operators 
are linear and may therefore be solved implicitly with ease. (The implicit- 


ness is necessary to overcome the stiffness which was transferred into the 
linear part.) The total scheme may be shown to be stable under the less re- 
strictive time step of the nonlinear part. In a future paper, we intend to 
present computational results for our proposed algorithm. 


1. A MODEL PROBLEM 

Consider the initial value problem for the following symmetric hyperbolic 
system 

“t-( u v), -(i ?) G) «•» 

' ' t X 

3 is a real number, |g| > 1. The eigenvalues of A are 


Uj(A) = 1 + e, p 2 (A) =1-3, 


and therefore an explicit scheme will have the CFL condition 
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For example, the Lax-Wendrof f scheme 
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is stable under the condition 


+ l*D <> 


Suppose now that one attempts to advance the solution of (1.1), equation by 
equation, rather than to use the form of the system as in (1.3). This amounts 
to splitting the matrix A into the sum of two matrices B and C 


A 





B + C 


( 1 . 4 ) 


and advancing the solution by using first the equation 


and then 
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where the initial value of (1.6) at every time step is the value of w^^ 
obtained after advancing (1.5) one time step. This procedure yields a scheme 
which is first order in time and second order in space. We note that the sys- 
tems defined in (1.5) and (1.6) are strictly hyperbolic and hence well- 
posed. The eigenvalues of B and C are 0 and 1, and therefore the Lax- 
Wendroff scheme for (1.5) and (1.6) separately will be stable under the condi- 
tion 

< 1 (1.7) 
Ax — 


allowing a time step much larger than the one allowed by (1.2) if 0 is a 
large number. However, even if a numerical method is stable for (1.5) and 

(1.6) separately, it need not be stable for the combination of (1.5) and 

(1.6) . In fact, consider the Lax-Wendroff scheme for (1.5) and (1.6). The 

amplification matrix G of the combined scheme is given by 
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where 
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kAx 
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We will show that the eigenvalues of G are greater than one in modulus for 
any A , and thus the combined scheme is unconditionally unstable. To do 
that we look at the mode £ = 1 : 


G(5 « 1) - 


n - 2\ + i 
' - 2 gx 2 


2 2 4 

■ 48 A 


-6(1 - 2A 2 )2A 2 


1 - 2A 


The eigenvalues of G y± are given by 


U± = 1 - 2X + 23 X 


2* 4 /TT 

± / 3 X 


+ 1 - 2X 


23 X 


The scheme is clearly unstable for 


since in this case y + > 1 for 3 > 1 and y >1 for 6 < 1. It 
is also easily verified that y + > 1 for 6 > 1 for any X. Thus, the 
splitting (1.5) - (1.6) is the wrong way of splitting. 

Perhaps a deeper insight is obtained if we Fourier transform (1.1), 
(1.5), and (1.6). The solution operator for (1.1) in Fourier space is 
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E(w,At) = e 


iAait 


where w is the dual Fourier variable. 

The solution operator for the split scheme (1.5), (1.6) over one time 
step is 

S(»,4t) . e lft»4t e ia,A1:_ (1.8) 


For every fixed o> 


S(co,At) t/At -► E(t) . 


However, since = C, = B, an expansion of the right ^nd side of (1.8) 

shows that 

S(a) ,At) = [I + B(e iuAt - 1)][I + C(e 1 “ At - 1)]. 


If we put Atco = tt , we get 


S(o) ,A t) = (I - 2B) (I - 2C) 


/4g 2 - 1 

\ -26 


26' 

- 1 , 


and for any 1 6 1 > 1 S(w,At) has eigenvalues larger than 1. This illus- 
trates the instability. 


2. THE ISENTROPIC EULER EQUATIONS 

The isentropic Euler equations in one space dimension may be written as 


w — 
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u 

— — 

u 
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u 

p 

t 
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= A w , 


( 2 . 1 ) 
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where u is the velocity, p is the pressure, p is the density, and y 
is the adiabatic constant of the fluid. The normalized equation of state for 
the fluid is 

P = P Y . (2.2) 


The eigenvalues of the matrix A in (2.1) are u - c and u + c, where 
c = is the sound speed. Thus, if we were to solve (2.1) by an explicit 

difference scheme, we would have to impose a CFL condition of the form 


At , const 
Ax — | u I +c 


(2.3) 


We wish to study (2.1) in the low mach number regime so that p = p^, 
where p^ is the base flow density. We define 


e = 


p 0 - p 


Then \e\ « 1. Using (2.2) we conclude 


P - P 0 = “ Y£Pq + 0(€ 2 ), 


where pq is the base flow pressure. 

One possible splitting for (2.1) [7], is to write A as the sum of two 
matrices Aj and A 2 as follows 



'0 

1/p o 


U 

1/p - 1 /p 0 

A = - 
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We then advance the solution of (2.1) by first using the equation 


»«> - *,»<», 

t lx’ 


(2.5) 


and then the equation 


( 2 ) . ( 2 ) 

"t ‘ ^2 w x ' 


( 2 . 6 ) 


Since is a constant matrix, we could solve (2.5) analytically thus doing 
away with any CFL restriction. The eigenvalues of the matrix A£, however, 
are - u ± i/y c Q e + 0(€ 2 ). Thus, the splitting (2.5 - 2.6) is not a 
hyperbolic splitting. 

To examine the stability of the split scheme, we examine the Fourier 
transform of the solution operator, S(At), over one time step. The Fourier 
transform of S is S: 


„ iA_o)At iA.uiAt 

S(w,At) = e e 


Let a = Cq wAt, and (3 = /y Cq € wAt. 


After some computation, we obtain 
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To first order in £ , we may write A 2 as 
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-y p 0 £ 
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We then have 


iA_wAt . .. 

2 iutoAt 

e = e 


coshg 


-i/y c 0 p Q sinhe 


-isinhg //y c Q p 0 


coshg 


Hence, 


rcoshgcosa - 


S(co,At) = e 


iuu)A t 


s Inhg s ina 
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HL 

C 0 P 0 


(coshgsina + sinhgcosa)n 


lipQCQ(/y sinhgcosa - coshgsina) /y sinhgsina + coshgcosa 


The eigenvalues of S(co,At) are roots of the polynomial 


p(A) = A^ - (2coshgcosa + sinhgsina )A + 1 = 0. 

H 

By Miller's criterion [8] the roots of p(A) are inside the unit disc if 
and only if 

0 = (coshgcosa + — — sinhgsina I 1 • 

2/*y 

If we let a = Cq a>At = ir , then 0 > 1, whenever g > 0. Hence, at least 

a 

one of the eigenvalues of S(a),At) lies outside the unit disc in a neigh- 
borhood Of a = 7T • 

If we were to solve (2.5) - (2*6) using a pseudospectral difference 
scheme, we would have to impose a CFL restriction of 


i 
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7T< 1 

Ax — c 


0 


to ensure the stability of our split scheme. For if we were to use the 

ilex 

Fourier modes {e }, k = -N,»»*,N, as a basis for our numerical scheme 

the mesh width Ax for the grid of collocation points would be given by 

a - 

Ax ~ W+T) * 

Since our stability condition dictates that 

c^NAt < it , 


the CFL condition for our scheme assumes the form 


At < _L_ 
Ax ~ C 0 * 


Nothing has been gained, therefore, from this splitting technique. 


3. THE EOLER EQUATIONS 


We write the Euler equations in one space dimension as 
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(3.1) 


Here p, m, p, and u denote the density, the momentum, the pressure, 
velocity, and c is the sound speed of the fluid. We analyze (3.1) in the 
low mach number regime. 

The eigenvalues of A are u - c, u and u + c. Therefore, an explicit 
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difference scheme for (3.1) would have a CFL restriction of the form 


At ^ const 


Ax — TuT+c 


(3.2) 


One possible splitting for (3.1) could be obtained by writing A as 



0 

1 
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O' 
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" A 1 + A 2 


(3.3) 


where Cq is the sound speed of the base flow. We then solve (3.1) by using 
first the equation 

W t° = A l w x° 0.4) 

and then 

f / o ^ 

(3.5) 


«i 2> - v< 2) - 


The advantage of such a splitting, it would seem, is that since A^ is a 
constant matrix we can obtain an analytical solution of (3.4) without any 
restriction on the time step. Further, since the eigenvalues of A2 are 0, 
-u, and -2u, we can solve (3.5) by a difference scheme with a large CFL con- 
dition of the form 

At / const f r\ 

T— " N VC! ”-T— ' . 

Ax — 2 |u| 


We examine the Fourier transform of the solution operator S(At) over one 

^ iA^rnAt iA.(*)At 

time step. Then S(w,At) = e 1 e 1 


Let 


a = c^wAt 
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2 2 
c- c Q 


We choose u = 0 and n > 0. Then 


lA.uiAt 

e =0 


isina 


-iCgSina 


1+cosa 

2 

C 0 
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iA-uiAt * ^ 

e = 0 1 

.1 -inc 0 a 
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Hence 


S(o),At) = 
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-i(ncQa cosa + c^sina) 


1+cosa 
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C 0 

isina 

C 0 
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The eigenvalues of S(o),At) are 1, and the roots of the polynomial 


p (X) = X -(2cosa - nasina) + 1 = 0. 


By Miller's criterion the roots of p(X) are inside the unit disc if and 


9 - 1 2cosa - nasina | ^ ^ 


only if 
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Let a = — 7T + 6 • 

2 

Then 0 = 1+ mr6 + 0(6 )• 

Hence S(w,At) has at least one root outside the unit disc near a = 

CgOdAt = it. Thus, this proposed splitting, once more, has undesirable proper- 
ties. If we were to solve (3.4) - (3.5) using a pseudospectral difference 

scheme, we would have to impose a very restrictive CFL condition of the form 

At . const 
Ax - c 0 

Gustafsson and Guerra [3] showed how to split (2.1) in a way that avoids 
the pitfalls pointed out above. The main idea in their work was to obtain two 
symmetric split operators. This, of course, is harder to do for more compli- 
cated systems. In the following sections, we generalize this approach to the 
problem of obtaining split operators which are simultaneously symmetrizable in 
the case of the full Euler and Navier-Stokes equations. 


4. CORRECT SPLITTING FOR THE EULER EQUATIONS 

In the preceding sections, we gave examples of "natural" splitting proce- 
dures which led either to instabilities or to stability conditions which at 
best did not represent an improvement over the original ones. A common fea- 
ture of those split operators was that they were not simultaneously symme- 
trizable. 

To avoid the dangers pointed out by these examples, we propose to remove 
the stiffness of a given stable symmetric operator by instituting a splitting 
procedure such that all the split-off operators are simultaneously symme- 
trizable. If each of these new operators is discretized in a stable manner, 
then the overall scheme will remain stable. 
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A prescription for a general operator achieving this goal is not known to 
us. We would like, however, to suggest such a procedure for compressible, low 
Mach number flows governed by either the Euler or the Navler-Stokes equa- 
tions. These systems are chosen in view of the "counter-examples" given in 
Section 3. The Euler equations may be symmetrized nonlinearly by using 

"entropy-variables" [5, 6J. The system thus obtained is of the form 

3 

p q t + I M* = 0 

i=l i 


where P and the A^s are symmetric matrix functions of the vector q. The 
premultiplying matrix P is usually non-sparse, and hence it is not clear how 
to remove the stiffness (if there is any) from the Aj,'s. In the Euler equa- 
tions, it is well known that the eigenvalues of are u, u + c, u - c 
where u is the x-component of velocity and c is the speed of sound. At 
low Mach number flows, u << c everywhere; hence, a Von-Neumann like stability 
condition 


At < 


Ax 

|u>c 


(4.2) 


gives an over-restricted condition. In this sense, the system is stiff (see 
Sections 2 and 3). 

Our approach is motivated by previous results [1] valid for the 
linearized frozen coefficient case. 

Consider the Euler equations for a gas in their nonconservative form in 
two-space dimensions (the three-dimensional case follows directly from the 
results of this section): 
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9t 9x 9y 


(4.3) 


A A A A A 

where V is the column vector (p,u,v,p) and the coefficient matrices 
are given by 


, (4,4) 


where p,u,v,p and y are, respectively, the density, the velocity com- 


ponents in the x and 


directions, the pressure and the ratio of 


specific heats at constant pressure and volume (y = C^/C^). Next, 
nondimensionalize the quantities in (4.3) as follows: 


T" ,x = T’ y = T ,p= F 


(4.6) 


where the subscript 00 indicates free stream condition and 


is a 


reference length. Equations (4.3) and (4,4) then retain the same form exactly 


with the superscript 


removed. In particular, the dimensionless speed of 


sound retains its functional form, i,e., 


c = — = /yp/p • 

u 


(4.7) 


We now propose the following change of variables: 
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-/y(y'-i) - 


(4.8) 


where is a constant to be specified later. One may then cast (4.3) in 

the form of (4.1) where: 
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(4.10) 


With these definitions of P, Aj , and A 2 the Euler equations 

Pq t + AjC, + A 2 q y - 0 


(4.11) 
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are symmetric hyperbolic* 

£ 

We now wish to motivate the way in which the operators in (4*11) (A^ 

£ 

and A 0 -r— ) are split. The starting point is the fact that we are here 
Z a y 

interested in low Mach number flow. Such flows are characterized by two 

2 2 2 

facts: the first is that c » u + v everywhere; secondly, for reason- 


able initial conditions 


2 , •_ \ 2 

c (x,y,t) - c 


« 1 . 


(4.12) 


For example, at steady state 


c (x,y) - c* __ T st - T. _ T .J u2 

o ' m o 


- T 


2 - 


(4.13) 


where T t is the stagnation temperature, M^ is the free stream Mach 

number; hence, for low Mach numbers (4.12) is valid. In view of the above, we 

choose c, = c , and we rewrite (4.11) as follows 
1 00 


Pq t + ^ R 1 + S P q x + ^ R 2 + S 2^ q y = ° 


(4.14) 


where 
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The four eigenvalues of P“^R^ are 

X = u, u, u ± (c - c ) [ 1 + — (— ) + — (— ) 2 ]^ 2 . 

/7 77 

It is clear from (4.11) that 


f 
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while 


c - C<a < (y-DM. 

C co 

— =■ 0 ( 1 ). 

C 


Thus, none of the eigenvalues gets to be large unlike the original unsplit 

scheme which had eigenvalues u, u, u ± c (recall that in our case u * 0(1) 

while c “ 1/M .) 

00 

Next we notice that and S 2 are constant matrices. A difficulty 

remains however in the nonlinear element of P. We shall deal with this as 
the method of solution is presented. 

Step I in the solution of (4.12) is to numerically advance 



The initial conditions for (4.16) are given by the solution to (4.15) at 
t = At. Notice that while S. and S ? are constant matrices, P has the 
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2 2 

nonlinear element c /c m . This means that removal of the stricter time-step 
due to Sj and S 2 (At _< const. Ax/cq) cannot be done easily via implicit 
method implementation of (4.16). To overcome this difficulty, we split 
and S 2 as follows: 


where 
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and 
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(4.18) 


Thus, we replace (4.16) with the sequence 


Pq. + S*q + sjq = 0 
^t lx 2 y 


(4.19) 


Pq + sPq + S^q « 0. 
^t lx 2 y 


(4.20) 
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Note that in (4.19), c t is identically zero, and so over that time step we 
take c « c(x,y) from (4.15). The rest of (4.19) is therefore linear (be- 
cause = c^(x,y) is known) and can be solved implicitly with relative 

ease. Alternatively, the first three equations in (4.19) may be combined into 
a variable coefficient wave equation for p, namely: 

,2 4 , 

— =- (An p) = — 5 V (An p); (4.21) 

3t c (x,y) 

u and v are then obtained directly from the middle two equations of 
(4.19). In (4.20) it is £n p that does not change over the time step. 
The rest of the system is linear with constant coefficients and may also be 
cast into a wave equation for c: 


3 


2 

c 



(4.22) 


and again u and v are found directly from the middle two equations of 
(4.20). 

This completes the splitting method for the Euler equation. The temporal 
and spatial discretization depends on the particular problem. Straight 
splitting as described here will result in only first order accuracy in 
time. Alternating the order of solving between (4.15) * (4.19) + (4.20) 

to (4.20) * (4.19) + (4.15) will yield second order in time, see [2]. 
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5. EXTENSION TO THE NAVIER-STOKES EQUATIONS (N. S. CASE) 

The Navier-Stokes equations may be written as: 


Pq + A q + A„q = B.q + B.q + Cq + F + F. 
M t l M x 2 y Pxx 2 yy xy 1 2 


where q, P, Aj , and A£ are as defined in the preceding section, 
quantities on the right hand side are given by: 
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(5.5) 
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/Y(Y-l) 
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where the dimensionless viscous production function $ is given by 


* = “ 7^ u x + V y )2 + 2 * u x + V y ^ + ^ u y + v x ]2 ' 


(5.1) 
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(5.2) 

(5.3) 

(5.4) 
(5.6) 


(5.7) 
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We can now describe the solution method: after obtaining the "hyperbolic" 

solution (see equations (4.15) to (4.22)), we go through the following steps: 


1) From (5.2) to (5.6), it follows that q = 0, i.e., during the "viscous" 

A t 

integration p = p(x,y). 

2) multiply equation (5.1) by the matrix 


0 

0 

0 

.0 


0 

1 

0 

0 


0 

0 

1 

0 


0 

0 

0 

0 


The resulting two "viscosity split" equations for u and v are 


u = 


1 [(iu + u ) + 4 u ] 


t ~ Ri- lv T 


XX 


’ -1J- U. j 

yy j xy 


(5.8) 


V — y [V V + y V ]. 

t Re xx 3 xy 3 xy 


(5.9) 


These may be easily solved implicitly since they are linear p.d.e.'s with con- 
stant coefficients. 


(3) The last step is to solve the viscous part of the energy equation which 
may be cast in the form: 


3 

7t 


(c 2 ) 


Y 

2PrRep 


[V 2 (c 2 )] + 


1 * 

TETp * 


(5.10) 
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Note that $/p is a function of the squares of u x , u y , v x , v y , and 
p(x,y) only and may therefore be taken as known from the previous step. 
Equation (5.10) is then a scalar linear inhomogeneous heat equation which 
again may be easily solved by implicit methods. 

Note that all the operators in (5.8) and (5.9) may be taken to be 
stable (i.e., !•! < 1) in L^. In addition, because c^ > 0, F 2 > 0 

(5.10) is also stabilizable under the Lj norm for c^; this assures the 

L 2 stability for q^. 

Notice the total algorithm (4.15) (4.22) (5.8) -*■ (5.10) may be 

run partly in parallel thus enhancing its efficiency beyond the removal of the 
stiffness. Schematically, the tree of calculation may be shown as follows: 



Thus, if parallel processors are available, we run only three computations 


instead of five 
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